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We examine the interaction of two-dimensional solitary pulses on falling liquid films. We make use of the 
second-order model derived by Ruyer-Quil and Manneville [Eur. Phys. J. B 6, 277 (1998); Eur. Phys. J. 
B 15, 357 (2000); Phys. Fluids 14, 170 (2002)] by combining the long-wave approximation with a weighted 
residuals technique. The model includes (second-order) viscous dispersion effects which originate from the 
streamwise momentum equation and tangential stress balance. These effects play a dispersive role that 
primarily influences the shape of the capillary ripples in front of the solitary pulses. We show that different 
physical parameters, such as surface tension and viscosity, play a crucial role in the interaction between solitary 
pulses giving rise eventually to the formation of bound states consisting of two or more pulses separated by 
well-defined distances and travelling at the same velocity. By developing a rigorous coherent-structure theory, 
we are able to theoretically predict the pulse-separation distances for which bound states are formed. Viscous 
dispersion affects the distances at which bound states are observed. We show that the theory is in very good 
agreement with computations of the second-order model. We also demonstrate that the presence of bound 
states allows the film free surface to reach a self-organized state that can be statistically described in terms 
of a gas of solitary waves separated by a typical mean distance and characterized by a typical density. 



I. INTRODUCTION 

The dynamics of a liquid film falling down a verti- 
cal wall has been the subject of numerous studies since 
the pioneering works by the father-son Kapitza teamir— . 
Wave evolution on a falling film is a classical long-wave 
hydrodynamic instability with a rich variety of spatial 
and temporal structures and a rich spectrum of wave 
forms and wave transitions, starting from nearly har- 
monic waves at the inlet to complex spatio-temporal 
highly nonlinear wave patterns downstream. Such pat- 
terns are known to profoundly affect the heat and mass 
transfer of multi-phase industrial units. Reviews of the 
dynamics of a falling film are given in Refs. i-d 

For small-to-moderate values of the Reynolds number 
(defined typically as the ratio of the inlet flow rate over 
the kinematic viscosity), the falling-film surface is pri- 
marily dominated by streamwise fluctuations and can 
be considered as free of spanwise modulations, i.e. two 
dimensional^. As it has been observed in many ex- 
perimental studies^— and theoretical worksi^r— , under 
these conditions, the film free surface appears to be ran- 
domly covered by localized coherent structures, each of 
which resembling (infinite-domain) solitary pulses. These 
pulses are a consequence of a secondary modulation in- 
stability of the primary wave field. They consist of a 
nonlinear hump preceded by small capillary oscillations 
and can even appear at small Reynolds numbers (but 
above the critical Reynolds number for the instability 
onset which vanishes exactly for a vertical film). 

Solitary pulses are stable and robust and continuously 
interact with each other as quasi-particles through at- 
tractions and repulsions giving rise to the formation of 



bound states of two or more pulses travelling at the same 
speed and separated by well-defined distances. Bound- 
state formation of two or more pulses has been recently 
observed experimentally in the problem of a viscous fluid 
coating a vertical fiber—. In this case, the initial growth 
of the disturbances is driven by the Rayleigh-Plateau 
instability and inertia. Eventually, the interface is domi- 
nated by drop-like solitary pulses which continuously in- 
teract with each other and can form bound states. The 
interplay between, on the one hand the Rayleigh-Plateau 
instability and inertia which enhance the front capillary 
ripples, and on the other hand surface tension and vis- 
cous friction that suppress them, has a crucial effect on 
the pulse interaction dynamics and, thus in turn, on the 
distances at which bound states are formed. 

In the present study we examine both analytically and 
numerically the bound-state formation phenomena in a 
vertically falling liquid film. We use a two-field model 
for the local flow rate and the liquid free surface that 
contains terms up to second order in the long-wave ex- 
pansion parameter— ~—. Hence, this model includes the 
second-order viscous effects originating from the stream- 
wise momentum equation (streamwise viscous diffusion) 
and tangential stress balance, i.e. second-order contribu- 
tions to the tangential stress at the free surface. These 
terms, which have been ignored in all previous pulse in- 
teraction theories for film fiows^i, have a dispersive ef- 
fect on the speed of the linear waves (they introduce a 
wavenumber dependence on the speed) and they affect 
the shape of the capillary ripples in front of a solitary 
hump. More specifically, increasing the strength of the 
viscous dispersive effects leads to decreasing the ampli- 
tude of the capillary waves ahead of the hump. This 
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effect is amplified as the Reynolds number is increased 
and hence should play an important role on the selec- 
tion process that brings the pulses to be separated by 
well-defined distances. 

We start by carefully developing a rigorous coherent- 
structure theory for the second-order two-field model. 
The aim is to obtain a dynamical system describing the 
location of each pulse by assuming weak interaction be- 
tween pulses (i.e. the pulses are sufficiently far from each 
other and they interact through their tails only). Al- 
though the basic ansatz we use at the outset, i.e. a su- 
perposition of N pulses plus an overlap function, is a 
standard assumption in weak interaction theories (and 
originates from condensed matter physics where it has 
been used to describe particle-particle interaction), the 
way we implement it into the two-field model is highly 
non trivial. For example, the spectral analysis of the 
resulting linearized non self-adjoint operator describing 
soliton interaction requires a careful and rigorous study 
that has not been performed before for a two-field sys- 
tem. For a single equation, the generalized Kuramoto- 
Sivashinsky (gKS) equation, a careful and rigorous study 
was performed recently in Refs. HH, [2(], and Hl|. Here 
we appropriately extend this study to the second-order 
two-field model. 

Of course, coherent-structure theories have been for- 
mulated for many different systems in recent years (see 
e.g. for the review of some of the methodologies for 
the gKS equation). Even rigorous justification of the 
ordinary-differential equations describing the dynamics 
of the pulses has been provided recently for equations 
with a stable primary pulse (e.g. Refs. [23 and [241). How- 
ever, for the present problem, as well as for the gKS 
equation analyzed recently in Refs. [lfl [2(| and [U the 
pulses are inherently unstable with the zero eigenvalue 
of the linearized interaction operator embedded into the 
essential spectrum. This in turn suggests that the usual 
projection procedure used in previous studies, such as 
those on weak-interaction approaches for the gKS equa- 
tion (e.g. Refs.[5l.[25l andUg), cannot be rigourously jus- 
tified. Moreover, previous studies appear to be either 
incomplete or at times overlook important details and 
subtleties. For example, the structure of the spectra of 
the adjoint operator of the linearized equation for the 
overlap function in the vicinity of a pulse, a crucial step 
for performing projections, has not been analyzed care- 
fully. This is done here by considering the linearized 
interaction operator on a finite domain and by imposing 
periodic boundary conditions. In this way we are able to 
recover the eigenvalues and the corresponding eigenfunc- 
tions on an infinite interval in the limit of the periodicity 
interval tending to infinity. We then show that the null 
adjoint cigcnfunction has a jump at infinity, which in 
turn implies that the localized function in the null space 
of the adjoint operator given in Ref.H (Fig. 9.1(c)) and 
also postulated in Ref. [26| is erroneous (these misconcep- 
tions are also discussed in the recent study by Tseluiko 
et al2L). 



As we shall also demonstrate here the projections for 
the derivation of the system governing the pulse dynam- 
ics can be made rigorous by the use of weighted functional 
spaces. We are then able to obtain rigorously a dynam- 
ical system describing the time evolution for the pulse 
locations. All possible distances at which bound states 
are formed are found by investigating the fixed points 
of this dynamical system. Detailed statistical analysis 
of time-dependent computations with the second-order 
model then shows that the separation distances between 
neighboring structures that the system selects are in ex- 
cellent agreement with those predicted by the coherent- 
structure theory. Moreover, the time-dependent compu- 
tations with the second-order model elucidate the influ- 
ence of viscous dispersion: increasing viscous dispersion 
allows the pulses to get closer to each other thus decreas- 
ing the separation distances between neighboring coher- 
ent structures. 

The use of weighted spaces makes our solitary pulses 
spectrally stable and hence allows also for the derivation 
of a novel and highly effective numerical scheme that can 
be used to accurately track the pulse dynamics for suf- 
ficiently long times. Furthermore, even though we pri- 
marily focus on weak interaction between solitary pulses, 
we also give numerical results on the strong interaction 
case, i.e. when the pulses are sufficiently close to each 
other, revealing a peculiar oscillatory behavior. Although 
an appropriate strong interaction theory describing this 
phenomenon is still lacking, we show that considering or 
not the second-order viscous dispersion effects becomes 
crucial for the description of the strong interaction dy- 
namics. 

In Sec. [IT] we present the second-order model that in- 
cludes the viscous-dispersion effects. In Sec. [Hi] we de- 
velop a coherent-structure theory for the interaction of 
solitary pulses of the second-order model. In Sec. HVl 
we contrast our theoretical predictions for formation of 
bound states with time-dependent computations of the 
full system. Finally, we conclude in Sec. [V] 



II. SECOND-ORDER MODEL 

Figure [T] shows the problem definition. We consider a 
thin liquid film flowing under the action of gravity down 
a vertical planar substrate. The liquid has density p, 
kinematic viscosity v and surface tension a. A Carte- 
sian coordinate system (x, y) is introduced so that x is 
in the direction parallel to the wall and y is the outward- 
pointing coordinate normal to the wall. The wall is then 
located at y = and the free surface at y = h(x,t). 
The governing equations are the mass conservation and 
Navier-Stokes equations along with the no-slip and no- 
penetration conditions on the wall and the kinematic and 
tangential and normal stress balance conditions on the 
free surface. 

By introducing a formal parameter e representing a 
typical slope of the film, e ~ \d x h\/h, we can perform 
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FIG. 1. Sketch of the profile geometry for a two-dimensional 
liquid film flowing under the action of gravity down a vertical 
wall. 



a long-wave expansion of the equations of motion and 
associated wall and free-surface boundary conditions for 
e <C 1. This is based on the observation that because sur- 
face tension is generally large, the interfacial waves are 
typically long compared to the film thickness. This so- 
called "long- wave approximation" has been central to all 
thin-film studies (see e.g. Refs. H and H3) and the small 
parameter e is frequently referred to as the "long-wave" 
or "film parameter" . The long- wave approximation leads 
to a hierarchy of model equations, starting from a single 
highly nonlinear partial differential equation for the film 
thickness h, the so-called "Benney equation"— whose re- 
gion of validity is restricted to near-critical conditions, 
to equations of the boundary-layer type which are valid 
away from criticality^— . By neglecting terms of 0(e 3 ) 
and higher, the so-called "second-order boundary-layer 
equations" read: 

u x + v y = 0, (la) 
Re(u t + uu x + vu y ) = 3(1 + Weh xxx ) 

-\-Uyy — 2ll XX + [lia;|/lja:) 

(lb) 

u\ Q = v\ o = 0, (lc) 

Uy\h = 4hx(v>x\h) — v x \h, (Id) 
h t + q x = 0, (le) 

where x and y have been non-dimensionalized with the 
Nusselt flat- film thickness h jv, the streamwise and cross- 
stream velocity components, u and v, respectively, are 
both non-dimensionalized by the average Nusselt flat-film 
velocity, ujy = gh 2 N /2>u, where g is the gravitational ac- 
celeration, and time t is non-dimensionalized by h^/iiN- 
Here, q = J Q u dy is the streamwise flow rate while the 
notations (-)|o and (-)\h indicate that the corresponding 
function is evaluated at y = and y = h(x,t), respec- 
tively. The two dimensionless groups appearing in (|lb[) 
are the Reynolds number, measuring the relative impor- 



tance of inertia to viscous forces and the Weber number, 
measuring the relative importance of surface tension to 
gravity, defined as: 
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We note that elimination of the pressure from the cross- 
stream momentum equation constitutes the main ele- 
ment of the boundary-layer approximation: by neglecting 
the inertia terms in this equation, the resulting equation 
can be integrated across the film to yield the pressure dis- 
tribution in the film which in turn is substituted into the 
streamwise momentum equation. It is also important to 
note that the second-order contributions in the long- wave 
expansion are given by the two last terms in the right- 
hand sides of (fib)) and (fTd|^£. Hence, the second-order 
boundary layer equations include streamwise viscous- 
diffusion effects and second-order contributions to the 
tangential stress at the free surface. These terms play 
a dispersive role, i.e. they introduce a wavenumber de- 
pendence on the speed of the linear waves^. 

By combining the long- wave expansion with a weighted 
residuals technique based on Galerkin projection in which 
the velocity field is expanded onto a basis with polyno- 
mial test functions, Ruyer-Quil and Mannevillei£~— ob- 
tained the following second-order two-field model for the 
local film thickness and flow rate: 
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where lengths, time and velocities in ([1} have been 
rescaled using the following scaling due to Shkadov^ 



(x,y) i-> (KX,y), 
(u, v) i-> (u, k~ 1 v), 

t H> Kt. 



where K = We 1 ^ 3 and 
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(5) 
(6) 
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corresponding to a reduced Reynolds and viscous- 
dispersion number, respectively (note that all second- 
order/viscous-dispersion terms are gathered in the sec- 
ond line of (pajl ). For rj = we obtain the first-order 
model 18 , which, much like the second-order one, can be 
derived from the first-order boundary-layer equations us- 
ing the long-wave expansion and a weighted residuals 
technique. It should also be noted that the second-order 
model contains the same number of parameters as the 
second-order boundary-layer equations in ([TJ and hence 
this model retains the "structure" of these equations (in 
contrast e.g. with the boundary-layer theory of aerody- 
namics where the Reynolds number can be scaled away 
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from the boundary layer as the corresponding equations 
are "simpler" compared to full Navier-Stokes) . 

This first-order model is similar to the first-order 
model obtained by Shkadov 30 by combining the long- 
wave approximation and averaging the basic equations 
across the fluid layer (effectively, a weighted residual 
technique with weight function equal to unity) but with 
different coefficients. As was pointed by Ruyer-Quil 
and Mannevillei^ the second-order model in (J3)) corrects 
the shortcomings of the first-order model obtained by 
Shkadov^, the principal one being erroneous prediction 
of the instability onset, i.e. of critical and neutral quan- 
tities: (a) it yields a 20% error of the critical Reynolds 
number for an inclined film (once again, for a vertical 
plane the critical Reynolds number vanishes) and (b) 
the neutral curve (wavenumber for the instability as a 
function of the Reynolds number) is not in agreement 
with Orr-Sommerfeld; for this purpose the second-order 
viscous-dispersion terms are crucial. This is a point 
that was analyzed carefully by Ruyer-Quil and Man- 
nevilleii. These authors contrasted the neutral curve ob- 
tained from both the first- and second-order model to the 
one obtained from the Orr-Sommerfeld eigenvalue prob- 
lem. The comparison shows clearly that the second-order 
model is in better agreement with Orr-Sommerfeld than 
the first-order one. 

Finally, it is noteworthy that the parameters 8 and 77 
can be expressed as 

s = r- 1 / 3 (3Re) n / 9 , ^r- 2 / 3 ^) 4 / 9 , (8) 

where T = er/^pp 1 / 3 !/ 4 / 3 ) is the Kapitza number that de- 
pends on the physical properties of the liquid only. It can 
then be readily seen that the second-order contributions, 
controlled by 77, are expected to be more relevant as Re 
increases and for liquids with either small surface tension 
or large viscosity. 



A. Steady-state solutions: solitary pulses 

Solitary pulses are traveling-wave solutions propagat- 
ing at constant speed Co = co(S, 77), hence stationary in 
a frame moving with speed Co, and sufficiently localized 
in space. Introducing then in ^ the moving coordinate 
x — > x — cot and requiring that waves are stationary in 
this moving frame, yields: 
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where qo(x) and ho(x) are the stationary solutions for 
the local flow rate and free-surface shape, respectively. 
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FIG. 2. Stationary pulse profiles for different values of S and 
two different liquids: water (W) and water+glycerin (WG) 
for both the first-order (77 = 0) and the second-order (77 > 0) 
models. 



Integrating once (0?) yields q = cqIiq + a, where a is 
an integration constant that can be fixed by demanding 
that the free-surface height approaches the Nusselt flat 
film solution away from the solitary hump, i.e. = I as 
x — > ±00, giving 



qo(x) = c [h (x) - 1] + i, 



(10) 



which in turn is used to eliminate qo in (JSJx) . The result- 
ing equation is solved numerically with a pseudo-spectral 
scheme in a periodic domain [-L, L] in which we take the 
Fast Fourier Transform (FFT) of the equation to obtain a 
system of nonlinear algebraic equations for the unknown 
FFT components of ho and the pulse speed cq. This 
system is solved using Newton's method by choosing an 
appropriate initial guess. 

Figure [2] shows several examples of stationary profiles 
h (x). We have used the values 5 = 0.98, 1.28, and 
1.82, which correspond to water (W) at Re = 3, 3.75, 
and 5, respectively, and to a water and glycerin 50% 
by weight mixture (WG) at Re = 1.59, 2.0, and 2.66, 
respectively (both were used as working fluids in the 
falling film experiments carried out by Liu et al^, but 
for inclined films with small inclination angle). The kine- 
matic viscosities of W and WG are vyj = 10 -6 m 2 /s and 
v wg = ^w, respectively, the densities are pw = 1.0 
g/cm 3 and pwc = 1-13 g/cm 3 , respectively, the surface 
tensions are aw — 69 x 10~ 3 N/m and awe = 72 x 10~ 3 
N/m, respectively, and the resulting Kapitza numbers 
are = 3364 and Ty/G = 334, respectively. It is worth 
noting that as long as 5 is kept fixed, the first-order model 
(7/ = 0) cannot distinguish, for instance, between W at 
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Re = 3 and WG at Re = 1.59. Therefore, the results for 
the first-order model correspond to both physical situa- 
tions. 

For 5 = 0.98, 1.28, and 1.82, we obtain the 
following values of the pulse speed and the maxi- 
mum height, {co,h m ): (1.15,1.19), (1.34,1.43), and 
(2.44,2.88), respectively, from the first-order model, and 
(1.15,1,20), (1.36,1.46), and (2.43,2.86), respectively, 
from the second-order model using W, and (1.18,1.23), 
(1.42, 1.53), and (2.4, 2.83), respectively, from the second- 
order model using WG. As it has been noticed by Ruyer- 
Quil and Mannevillei^, the differences between the first- 
and the second-order models as far as solitary pulses are 
concerned become more significant as 5 is increased. Al- 
though the pulse speed and the maximum height do not 
differ much (differences are appreciable at low values of 
S), the capillary ripples that are present downstream of 
the hump appear to be largely suppressed as a conse- 
quence of the second-order viscous-dispersion effects, spe- 
cially for WG which has higher viscosity. 



III. COHERENT-STRUCTURE THEORY FOR 
INTERACTING PULSES 

As it has been emphasized in the Introduction, at suf- 
ficiently large distances from the inlet of the film, the 
dynamics of the free surface is dominated by the pres- 
ence of localized coherent structures, each of which re- 
sembling (infinite-domain) solitary pulses, which interact 
indefinitely with each other as quasi particles through at- 
tractions and repulsions. The objective of this section is 
to appropriately extend the recently developed coherent- 
structure interaction theory for the solitary pulses of the 
gKS equatio n 16 i 20 ' 21 to the two-field model given by 
This is by far a non-trivial task as we shall demonstrate, 
e.g. unlike we now deal with a rather involved vector 
equation instead of the scalar gKS equation i n 16 i 20 ' 21 . 
The aim is to obtain a dynamical system describing the 
location of each pulse by assuming weak interaction be- 
tween pulses (i.e. the pulses are sufficiently far from 
each other and they interact through their tails only). 
This concept has been applied in many other fields, such 
as particle physics and quantum mechanics, where one 
typically deals with a system compound of many parti- 
cles, and has been used successfully to describe particle- 
particle interaction. As emphasized in the Introduction, 
previous coherent-structure theories appear to be either 
incomplete or some times overlook important details and 
subtleties, for instance in relation of the spectrum of 
the operator describing interaction. Moreover, previous 
coherent-structure theories for film flows in the region of 
small-to-moderate Reynolds numbers^^ were based on 
the Shkadov model and thus ignored the effect of viscous 
dispersion. 

We start by considering the system of partial differen- 
tial equations ([3]) in the frame moving with velocity cq of 
a single stationary pulse, x — > x — c$t. Then we assume 




FIG. 3. Schematic representation of a solution consisting of 
a supersposition of N pulses located at Xi for i = 1, . . . , N. 



that a solution for the local flow rate, q(x,t), and free- 
surface profile, h(x,t), is given as a superposition of N 
quasi-stationary pulses located at x\(i), . . . , ccjv(i) and 
a small overlap function, i.e. we postulate the following 
ansatz: 

1 N 

q(x,t) = -+y2Qi(x,t) + Q(x,t), (11a) 
6 i=i 

N 

h{x,t) = l + Y^ H i{x,t)+H{x,t), (lib) 

i=l 

where Qi(x, t) = qi(x, t) — 1/3 and Hi(x, t) = hi(x, t) — 1 
and 

qi(x,t) = q (x-Xi(t)} hi(x,t) = h (x-Xi(t)), (12) 

with qo and ho denoting the steady-state solution of 
© . A schematic representation of an iV-pulse solution is 
given in Fig. [3j where we have defined the distances be- 
tween two pulses as ii = x; l+ \ — Xi, for i = 1, . . . , N — 1. 
Next, we consider weak interaction between the pulses 
by assuming that they are well separated, i.e. li 3> 1, so 
that for each pulse it is enough to consider its interaction 
with only the immediate neighbours. We define a small 
parameter e — exp (— mini{£;}), and we assume that the 
pulse velocities x(t) and the overlap functions Q and H 
are of O(e). By substituting (fTTj) into (|3|) and expanding 
up to 0(e), we obtain a linearized equation in the vicinity 
of the ith pulse for the overlap vector function 



that is written as follows: 



(13) 



where 



n t -Xi(t)& ix =C i Sl + J i Y i , (14) 

Qi-l + Qi+1 



qi 

hi) 1 



Hi-i + H i+ i 



(15) 



for i = 1, . . . , N. Note that for i = 1 and i = N, 
the vector Y reads as: Ti = (Q,2,H2) T and = 
(Qn-i, Hn-i) t . d and J \ are the following linear ma- 
trix/differential operators: 



J- - \ J l J i) , (16) 
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Equation (fT4")) reveals that the dynamics of the overlap 
function in the vicinity of the ith pulse depends on the 
spectrum of the linear operator Ci and is influenced by 
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the neighbouring pulses, as indicated by the last term on 
the right-hand side of 



A. Analysis of the structure of spectrum of the linearized 
operator describing interaction 

It can be verified numerically that on a periodic do- 
main the operator Ci has a zero eigenvalue with geomet- 
ric multiplicity one and algebraic multiplicity two (the 
numerical scheme for constructing the spectrum of Ci 
will be described shortly). The operator d then has 
a null space spanned by the eigenfunction &\ = & ix 
that is associated with the translational invariance of the 
system. The corresponding generalized zero eigenfunc- 
tion 3> 2 that satisfies Ci&\ = &\ is associated with the 
Galilean invariance of the system. The aim then is to 
project the dynamics of the overlap function onto the 
null space of Ci in the vicinity of the ith pulse. 

From a physical point of view, the existence of these 
two dominant modes means that any perturbation to the 
steady pulse solution will make the pulse to shift (due to 
the translational mode) and/or to accelerate (due to the 
Galilean mode). We note that according to the solution 
ansatz (jlll) . we can assume that the overlap function, Si, 
is "free of translational modes" . The precise meaning of 
the latter phrase will be explained later. 

Projection onto the null space of the linear ma- 
trix/differential operator Ci requires a careful and de- 
tailed analysis of its spectrum as well as the spectrum of 
the adjoint operator C* (see Appendix |A")) . The essential 
spectrum of the operator Ci is given by the dispersion 
relation of the basic Nusselt state, (ftv,hjv) = (1/3,1), 
in the moving frame. By replacing (ft, hi) in Ci with the 
Nusselt state, Ci becomes a matrix operator with con- 
stant coefficients and its essential spectrum X(k) satisfies 
the following equation: 



+ ( c o 



§(i£0 2 -A A + i iA:+ ^ (ifc) 3_| i(ifc) 
coifc — A 



(17) 



for fcel. The locus of the essential spectrum is shown 
as a solid line in Fig. |4l and coincides with the locus of 
the essential spectrum of the adjoint operator (see Ap- 
pendix [2} , as expected. We note that part of the essen- 
tial spectrum is unstable. As it has been demonstrated 
in previous studies (e.g. Refs. [13, US and l2ll). the un- 
stable part of the essential spectrum is connected with 
the flat film instability and can be excluded from our 
consideration. 

To analyze numerically the spectrum of Ci, we first 
compute the solution for hi and qi on a finite 2L-periodic 
interval by using a pseudo-spectral method. The ma- 



trix representation of the linear operator is obtained by 
applying each component of Ci to a plane wave, i.e. 
L\ — > £je lk ™ x , for j = 1, 4, and transforming to 
Fourier space the resulting functions. By computing this 
operation Vfc„ = nir/L of the truncated Fourier series 
of hi, we are able to get the nth column of the Fourier 
matrix representation of C\ . We then analyze the eigen- 
values and eigenvectors for the resulting matrix. Our 
results show that, in addition to the essential spectrum 
and the zero eigenvalue, which is embedded into the es- 
sential spectrum, there is one more eigenvalue, which is 
negative and isolated. The eigenvalues are depicted as 
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FIG. 4. Spectrum of d for 8 — 0.98 for the second-order 
model with the physical parameters corresponding to W. The 
solid line represents the essential spectrum and the crosses 
represent the point spectrum. 
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FIG. 5. Components of the generalized zero eigenfunction 
<I?2 = (09i)<^hi) corresponding to the Galilean mode. 



crosses in Fig.@] In Fig. [S] we depict the two components 
of the generalized eigenfunction <& l 2 given by the solution 
of Ci&2 = &\ and corresponds to the Galilean mode. 

On a periodic domain, the discrete part of the adjoint 
operator also coincides with the discrete part of Ci (see 
AppendixEJ ■ We find that the eigenfunction correspond- 
ing to the zero eigenvalue on a periodic domain is merely 
a constant 



(18) 



and the generalized zero eigenfunction *$> 2 = (ip qi ,if)hi) T 
has to be found numerically. Its components are shown 
in Fig. [6] for different values of the periodicity interval. 

In our analysis, we have normalized the eigenfunctions 
so that the following conditions hold: 



<*i,*i> 



= o, 

= i, 



(* 2 ,* 2 )=0, 



where (-, •) denotes the inner product in L^(—L, L). 

We now examine the behavior of the eigenfunctions in 
the limit L oo. First, we note that both components 
<f> qi and <phi corresponding to the Galilean mode <fr 2 do 
not decay to zero at infinities x — > ±oo (cf. Fig. [SJ, and 

therefore we have that \\& 2 \ = J(& 2 ,& 2 ) -> oo as the 
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FIG. 6. Components of the generalized adjoint eigenfunction 
^2 = i^qn^hi) for different values of the period L. 




FIG. 7. (a) Norm of the generalized zero eigenfunction of 
d, for different values of the system size L. The solid line 
corresponds to a function proportional to L 1 ^ 2 . (b) Constant 
value m of the zero eigenfunction of £* for different L. The 
solid line corresponds to a function proportional to 1/L. 



system size is increased, see Fig.JT^a). This means that on 
an infinite domain the zero eigenvalue has both algebraic 
and geometric multiplicity equal to unity and the null 
space of Ci is spanned only by the translational mode. 
Second, we also note that the constant m in (fT5|) corre- 
sponding to the zero eigenfunction of C* tends to zero as 
the system size is increased (see Fig.[7][b)), meaning that 
C*^f l 2 — > as L — > oo, and thus, we have that the func- 
tion ty 2 belongs to the null space of C*. As it is shown 
in Fig. [SI the component i\) qi decays to zero at infinities 
for L — > oo, and the component ipht tends to different 
constants as x — > ±oo, and therefore \\*& 2 \\ oo. 

We then conclude that the zero eigenvalue of Ci is 
not isolated but belongs to the essential spectrum (and 
has both algebraic and geometric multiplicity equal to 
unity on an infinite domain with the null space of Ci 
spanned by the translational mode <&\). Also, zero is not 
in the point spectrum of the adjoint operator, C* on an 
infinite domain (its null space is spanned by a constant 
and a vector function, one component of which does not 
decay to zero at infinities). These two points complicate 
the formal projections of the overlap function onto the 
translational mode (the null space of Ci). 

We remark here that we find similar spectrum fea- 
tures between the second-order model and the gKS 
equatio n 16 i 20 ' 21 . We first note that although the Galilean 
mode $2 m the gKS equation is given by a constant, it 
is observed in both models that its norm tends to infin- 
ity as L is increased, and therefore it can be excluded 
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from the projections on the translational mode. In ad- 
dition, we find the same behaviour of the /i-component 
of the generalized eigenfunction of the adjoint operator 
(cf. Fig. and the generalized eigenfunction of the in- 
teraction problem with the gKS equation, namely a jump 
at infinity in both cases and, therefore, an infinite norm 
for the corresponding eigenfunction. As was pointed out 
in the coherent-structure theory for the gKS prototype^, 
it is possible to overcome such difficulties by making use 
of a formal procedure in a weighted space of functions 
that decay exponentially at +00. We shall therefore ap- 
ply such a formalism in the present study. 

B. Formulation in a weighted functional space 
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FIG. 8. Essential spectrum of CI for S = 0.98 with physical 
parameters for W and different values of a. 



Projections onto the null space can be made rigorous 
by choosing an appropriate weighted space. Following 
the formulation used in |3ll for the Korteweg-de Vries 
equation and in [2(3 andl2ll for the gKS equation, we shall 
restrict our projections to the following weighted space: 



Ll = {u : e°*u G L 2 }, 



(19) 



with the inner product (u, v) a — (e ax u, e ax v), and a > 0. 
The spectrum of the matrix/differential operator Ci in 
L\ can be studied there through the matrix/differential 
operator defined by 



Clu = e ax Ci{e- ax u), 



(20) 



on Lg. More precisely, one can easily see that the es- 
sential spectrum X a (k) of £" is then given by ([T71) by 
replacing ik with ifc — a. The interesting point of working 
in such a weighted space is that for certain values of a, it 
is possible to completely shift the essential spectrum to 
the left in the complex plane (see Fig. [8]), and therefore, 
the pulses may become spectrally stable by choosing an 
appropriate value of a. Such a shift of the essential spec- 
trum means that the pulses are transiently unstable but 
not absolutely unstable, i.e. any localized disturbance is 
convected to — oo in a frame moving with the velocity 
of the pulsed. [The essential part of the spectrum typi- 
cally leads to an expanding and growing radiation wave 
packet such that its left end travels upstream and right 
end downstream (absolute instability). If the speed of 
the right end of the wave packet is larger than that of a 
solitary pulse, the pulse is destroyed, i.e. it is "absolutely 
unstable"]. 

It is this absolute instability of the pulses which is re- 
sponsible for the complex turbulent-like spatio-temporal 
chaos observed in the KS equation, i.e. disordered dy- 
namics in both time and space and without clearly identi- 
fiable soliton-like coherent structures. For the gKS equa- 
tio n 20 i 21 there exists a critical value of the parameter con- 
trolling dispersion (the coefficient multiplying the third- 
derivative term in the equation), below which it is no 
longer possible to shift the spectrum to the left half-plane 
and the behavior of the gKS equation is spatio-temporal 



chaos like with the KS one. In our case, however, for all S 
values (or, equivalently, the Reynolds number) we exam- 
ined (up to 12) it is always possible to shift the spectrum 
completely to the left half of the complex plane, implying 
that dispersion effects prevent the system from evolving 
into spatio-temporal chaos. 

From the numerical point of view, having a stable es- 
sential spectrum means that the instabilities due to nu- 
merical noise could be eliminated and the temporal evo- 
lution of solitary waves could be obtained for sufficiently 
long times. In order to have then a stable essential spec- 
trum, we substitute 

q = e- ax g+l/3, h = e- ax f + l (21) 

into ([3]) to obtain the following equations for g and /: 

( ,5, 5 g - (l/3)e-°* f 2 - (2/3) / 
9t = c ( 9x - ag) + -f - 



-3af xx + 3a 2 f x - a 3 f) + 



^c- ax (f x -af) 2 



9 e~ ax q 
~2~~h~ (y9x ~ ag ^ x ~ a ^ ~ 6 fr(f xx 
g 

-2a f x + a 2 f) + - (g xx - 2ag x + a 2 g) 



ft = co(f x - af) - (g x - ag). 



(22a) 



(22b) 



Integrating numerically the above equations with an 
appropriate value of a, ensures that the flat-film insta- 
bility will not be excited and the dynamics of the pulses 
will not be affected by any numerical noise. This will 
be particularly useful, for instance, when studying the 
interaction between two pulses, see Sec. IIV A[ 

The other important consequence of working in a 
weighted space is that the projections onto the null space 
can now be made properly, since the zero eigenvalue be- 
comes isolated. We note that the null space of LI is only 
spanned by the zero eigenfunction, <fr l a = e ax &\. Also, 
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the zero eigenfunction of the adjoint operator, defined by 
C a* = e -^c*{e ax f), can be written as 

*i =e -«»(**- lim <J>|) 

x— > — oo 

which now has a finite norm. Therefore, we will use the 
projection operator 

Pi(f) = {f,K)a*l (23) 
for projecting onto the the null space of 

C. Pulse interactions and bound-state formations 

By applying the projection operator Pj to (] 14[) rewrit- 
ten in an exponentially weighted space and assuming 
that the overlap function is "free of translational modes" 
meaning that it is in the null spaces of the projections, 
i.e. Pi(e ax tt) — 0, we obtain the following dynamical 
system for the pulse locations: 

/oo 
[(#<_! + H i+1 ){c Jl* + J?*)*P gi ] ds, (24) 
-OO 

for i = 2, . . . , N — 1, and 

/oo 
[H 2 {c Jl* +J?*)ip qi ]dx, (25) 
-oo 

/oo 
[Hjr-^coJ}? + J^ qN ] dx, (26) 
-oo 

where Jl* and Jf* correspond to the adjoint operator 
components of the matrix/differential operator J7"i (see 
Appendix |A|) . and ip qi is the first component of the ad- 
joint eigenfunction Here we have also made use 
of the relation Qi — cqHi. If we define the function 
ipi = {cqJI* + Ji*)ipqii and use the notations 

/oo 
H (x-£/2)^ (x+e/2)dx, (27) 
-oo 

/oo 
H Q (x + £/2)^ (x-£/2)dx, (28) 
-oo 

where Hq(x) = ho(x) — 1 and ip is such that ipi( x ) — 
ip (x — Xi), ([24)) can be finally rewritten as 

±i(t) = Si(x i+ i - Xi) + S 2 (xi - Xi-i), (29a) 

for i = 2, N — 1. For i = 1 and i = N, we have 

ii(t)=5i(a:2-aji) (29b) 

and 

= ^(iCjv - Xff-i), (29c) 

respectively. Therefore, the time-evolution of the ith 
pulse location is controlled on the one hand by the 
function Si(x), which describes the interaction with the 



monotonic tail of the downstream pulse i + 1, and on 
the other hand, by the function S^x) that describes the 
interaction with the oscillatory capillary ripples of the 
upstream pulse i — 1. 

Let us consider the case of only two pulses interact- 
ing with each other, i.e. a binary interaction scenario. 
From (|29a|> . we can write an equation for the separation 
distance between the pulses, i(t) — x%{b) — x\(t), as: 

l(t) = S 2 (£)-S 1 (£). (30) 
The fixed points of the above equation are given by 

Si(tn) = S 2 (tn), (31) 

which predicts the different distances £ n at which both 
pulses travel at the same velocity, giving rise then to the 
formation of bound states. Figure |9] shows the graphs 
of Si and S 2 for two different values of S for the second- 
order model using W. Interestingly, since Si (£) and S 2 (£) 
represent the velocity of the pulses located at X\ and x 2 , 
respectively, we can also predict both the velocity of the 
bound state relative to Co, i.e. c n = Si(£ n ), and its sta- 
bility. Stable and unstable bound states are represented 
by solid circles and crosses, respectively in Fig. [9j As 
long as S 2 > Si, we have that x 2 > x\ so that the sec- 
ond pulse moves faster than the first one leading to an 
increase oil, and therefore, both pulses repel each other. 
On the other hand, when S 2 < Si, the first pulse is mov- 
ing faster than the second one leading to a decrease of £, 
and therefore both pulses attract each other. 

From a physical point of view, such oscillatory behav- 
ior of attractions and repulsions between pulses can be 
understood in terms of the interaction between the pe- 
riodic capillary waves ahead of the first pulse and the 
monotonic tail behind the second pulsed. More pre- 
cisely, the oscillatory shape of the free surface ahead of 
the hump induces a sign change of the capillary pressure 
difference across the free surface, Ap c ~ crh xx . When 
the tail of the second pulse overlaps with a maximum 
of one of the capillary waves of the first pulse, there is 
a drainage process of liquid from the first to the second 
pulse due to the positive pressure difference. This, in 
turn, generates an increase of both the height and the 
speed of the second pulse, on the one hand, and a de- 
crease of the height and the speed of the first pulse, on 
the other hand. As a result, the distance between the 
pulses increases corresponding to repulsion. The oppo- 
site behavior occurs when the tail of the second pulse 
overlaps with a minimum of one of the capillary waves of 
the first pulse, giving rise then to an attraction process. 
When the frequency and the amplitude of the capillary 
oscillations are increased by increasing 8, the number of 
bound states observed in a given interval by such at- 
traction/repulsion mechanism increases accordingly [see 
Fig. Etb)], as expected. 
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FIG. 9. Functions Si and S2 (solid and dashed lines, respec- 
tively) for (a) 5 = 0.98 and (b) 5 = 1.82. The black circles 
and crosses correspnd to stable and unstable bound-state sep- 
aration distances, respectively. 



IV. NUMERICAL RESULTS 

In this section we present extensive numerical experi- 
ments for the second-order model ©■ More specifically, 
we investigate numerically the temporal evolution of a 
superposition of two pulses as an initial condition and 
resulting attractive and repulsive dynamics giving rise to 
formation of bound states. We contrast the numerical 
results with the coherent-structure theory developed in 
the previous section. By imposing a localized random 
initial condition, we are able to study how several pulses 
interact with each other to self-organize and form bound 
states compound of two or more pulses. The effect of 
viscous dispersion on the bound-state formation will be 
elucidated by systematically integrating both the first- 
and second-order models for the physical parameters cor- 
responding to W and WG. 



A. Superposition of two pulses 

Two-pulse bound states can be constructed numeri- 
cally by solving (jH]) with an initial guess consisting of 
a superposition of two pulses separated by the theoreti- 



cally predicted distance. We use the numerical method 
of Sec. Ill At based on a combination of a pseudo-spectral 
method and Newton's method, and the numerical result 
is compared to the superposition of two pulses, 



1 + H (x + in/2) +H (x- in/2), 



(32) 



where £ n corresponds to the distances obtained by pip . 
Figure [TU] shows the results in the case of 5 — 0.98 for the 
bound states predicted at £\ ss 13.9, £2 ~ 17.3, £3 20.7, 
and £4 « 27.5 [cf. Fig. Ufa)]. As expected, the longer 
separation distance between the pulses is, the better the 
agreement between the theory and the numerical solution 
becomes. In particular, at £3 20.7 and £4 « 27.5, the 
numerical solution and the superposition of two pulses 
with the theoretically predicted bound-state separation 
distance are practically indistinguishable. It is interest- 
ing to note that the numerical scheme only converged 
when the pulse separation distance for the initial guess 
was sufficiently close to the value predicted by the theory. 
We also find good agreement between the bound-state 
velocities relative to cq predicted theoretically and found 
numerically: the numerical values are c\ ~ —0.00302, 
c 2 fa -0.00039, and c 3 -0.00006, which are to be 
contrasted with the theoretical predictions of —0.00197, 
-0.00035, and -0.00006, respectively. 

To check the attraction/repulsion dynamics between 
two pulses predicted by the theory, we study numerically 
the time-evolution of two pulses separated by an initial 
distance £0 which is close to either a stable or unstable 
bound-state separation distance. We numerically inte- 
grate both the second-order model (|22l) in a weighted 
space, and the theoretical model given by (|30|) . As it has 
been emphasized in Sec. 1111 B[ working in the weighted 
space l? a has the advantage of the essential spectrum be- 
ing shifted to the left half of the complex plane. Hence, 
any instability on the flat film region that originates from 
numerical noise is eliminated, which then allows us to 
follow the temporal evolution of the two pulses for suf- 
ficiently long times. Equations (|2"2"|) are integrated by 
choosing a = 0.1. To solve (f2"2"j) numerically we use the 
FFT to obtain the Fourier and the inverse Fourier trans- 
forms of g and / in the right-hand sides of (|2"2"|) and a 
fourth-order Runge-Kutta method to march forward in 
time. A typical time step is At = 0.0025 and the periodic 
domain [—L,L], with L — 120, is discretized into 2000 
intervals. 

The results for S — 0.98 and 1.82 for W are presented in 
Fig- LIU The interaction between the two pulses is attrac- 
tive when the pulses are initially separated by £$ = 29.5 
and 36.4 for 5 — 0.98 and 1.82, respectively, converging 
to the predicted stable bound states with the separa- 
tion distances approximately 27.5 and 35.3, respectively. 
Likewise, the dynamics is repulsive when the pulses are 
initially separated by £0 — 33 and 36.9 for S = 0.98 and 
1.82, respectively. The separation distances then con- 
verge to the predicted stable bound-state separation dis- 
tances 34.3 and 37.5, respectively. In all cases, we found 
a very good agreement between the theory and the nu- 
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FIG. 11. Evolution of the separation distance between two pulses for 8 — 0.98 and 8 = 1.82. The initial separation distances 
are £o = 29.5 (a) and 33 (c) for 8 = 0.98, and lo = 36.4 (b) and 33.9 (d) for 8 = 1.82. Solid lines correspond to the numerical 
solution of the theoretical model given by (|30p . and dashed lines are the numerical solution of the second-order model in a 
weighted space, (|22[) . with a = 0.1. 



merical results. 

Interestingly, a different interaction behaviour emerges 
as both pulses are placed close enough to each other 
so that the assumption of a well-separated distance be- 
tween them is not valid anymore and the interaction is no 
longer weak. Figures [T27a) and [FZ[ c) depicts the tempo- 
ral evolution of two pulses separated by an initial dis- 
tance £q = 20 for a relatively high Reynolds number 
(5 = 1.82) in the second-order model using the physi- 
cal parameters for W (which gives a value of rj = 0.015 
for the viscous-dispersion number) . We observe that both 
pulses attract and repel each other in an oscillatory man- 
ner, giving rise to a rapid fluctuating growth of the initial 
separation length until they reach an oscillatory steady- 
state of constant frequency and amplitude [the dashed 



line in Fig. [T"2T c)]. Is is remarkable that although such 
an oscillatory dynamics cannot be captured by our weak- 
interaction theory, the final separation length around 
both pulses are oscillating corresponds to a stable bound 
state predicted by the theory (solid lines), and the am- 
plitude of the oscillations is delimited within the distance 
between two consecutive unstable bound states (dotted 
lines). This type of oscillatory interaction between two 
pulses was first observed by Malamataris et al>2 3 in di- 
rect numerical simulations of the full Navier-Stokes equa- 
tions with wall and free-surface boundary conditions, and 
was attributed to the competition of the strong fluctu- 
ations on the capillary pressure at the overlapping area 
between the pulses, and the nonlinear response of the 
solitary hump when it is perturbed from its stationary 
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FIG. 12. Oscillatory interaction between two pulses observed at S — 1.82 in the second-order model by using the physical 
parameters of W [(a), (c)] and WG [(b), (d)]. The top panels show the two-pulse system profile at four different times: t = 0, 
250, 400, and 600, from bottom to top, respectively. The bottom panels show the comparison between the numerical result 
obtained by integrating the second-order model in a weighted space, given by equations (1221) with a = 0.1 (red dashed line) 
and the asymptotic result predicted by the theoretical model (|30|) . The solid and dotted lines in (c) represent the stable and 
unstable bound states (BS), respectively, predicted by the theory. The pulses seem lock on at certain average distances close 
to those predicted by the weak interaction theory. 



shape. In this sense, if we consider stronger viscous dis- 
persion effects, leading to smaller amplitude of the capil- 
lary fluctuations in the front tail of the pulse (cf. Fig. [5]) , 
such oscillatory interaction can be largely or completely 
removed. Indeed, Fig. repeats the same numeri- 

cal experiment but taking the physical parameters corre- 
sponding to WG, which is more viscous than water and 
thus gives the value of rj = 0.052 for the viscous dis- 
persion number. We observe that such oscillatory inter- 
action is largely reduced and the two pulses eventually 
reach a stable bound state in agreement with our weak- 
interaction theory. This is a clear manifestation of the 
fact that including the second-order viscous-dispersion 
effects, ignored in all previous film- flow interaction stud- 
ies, can be crucial to obtaining an accurate description 
of the interaction between pulses. 



B. Localized random initial condition 

We also integrate ([3]) by considering both the first- 
(j] = 0) and second-order model (rj > 0) , and imposing a 
localized random initial condition. In our simulations, we 
have used the parameter values for both W and WG. To 



obtain an appropriate random representation of a differ- 
entiable function but with a sufficiently high frequency 
content we construct the initial condition by using the 
following random function: 



h in (x) 



To 



N 

E 

m — 1 



, sin ——X+Pr, 



kom 

cos x 

N 



(33) 



which has been recently proposed in |34| in the context of 
wetting of disordered substrates. Here, 70 and k are the 
characteristic amplitude and wavenumber, respectively, 
TV is a large positive integer, and a m and /3 m are sta- 
tistically independent normal variables with zero mean 
and unit variance. It can be shown that in the limit of 
7o^o 1 an d N — > 00, hi n (x) is a band-limited white 
noise. To get a localized perturbation, h- m (x) is then mul- 
tiplied by the function 6(x) = [tanh(a;) — tanh(a; + L c )]/2, 
where L c is the length of the disturbance. In our simu- 
lations we have chosen 70 = 0.05, fco = 1, N = 2000, and 
L c = 500. 

To accelerate the computational efficiency of our nu- 
merical scheme, we have used a Runge-Kutta-Fehlberg 
method with dynamic time-step adjustment. At each 
time step this scheme computes two different solutions 
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FIG. 13. Typical numerical solutions of the second-order model ((3| in extended domains for 5 = 0.98, 1.28, and 1.82 at 
t — 2250, 1700, and 450, respectively (panels (a), (b), and (c), respectively). 
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FIG. 14. Histograms of the pulse-separation distances for S — 0.98 in (a) the first-order model (rj — 0), (b) the second-order 
model with the parameter values corresponding to W (?) = 0.011), and (c) the second-order model with the parameter values 
corresponding to WG (r) — 0.041). The table shows the locations of the peaks observed in the histograms obtained for the 
lst-order (IstO) and 2nd-order (2ndO) models compared to the values predicted by (|29a[) . plotted in the insets of each panel, 
where the dashed and solid lines correspond to the S2 and Si functions, respectively. 
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FIG. 15. Histograms of the pulse-separation distances for S = 1.28 in (a) the first-order model (rj — 0), (b) the second-order 
model with the parameter values corresponding to W (rj = 0.013), and (c) the second-order model with the parameter values 
corresponding to WG (rj = 0.046). The table shows the locations of the peaks observed in the histograms compared to the 
values predicted by the model given by ([29a} , plotted in the insets of each panel, where the dashed and solid lines corresponds 
to the S2 and Si functions, respectively. 



and compares them. The time-step is the redefined ac- 
cording to the agreement between both solutions. We 
impose a minimum time step At = 0.001 (that is reached 
only in computations for large values of 6), and we are 
able to use time steps as large as At = 0.01. We inte- 
grated equations ([3]) on a periodic domain [-L, L] dis- 
cretized into 2M intervals. We used the following values: 
L = 1250, 1800, and 1200 and M = 2000, 3000, and 4000 
for 5 = 0.98, 1.28, and 1.82, respectively. The initial dis- 
turbance results in an expanding wave packet whose left 
envelope travels upstream (absolute instability; see our 
comment in Sec. IIIIBI) and the right one travels with a 
velocity lower to that of an individual pulse. The wave 
packet grows as it propagates and gives birth to a number 
of pulses escaping the expanding wave packet. The equa- 
tions were integrated up to t = 2250, 1700, and 450 for 
<5 = 0.98, 1.28, and 1.82, respectively. Beyond these times 
the rear side of the expanding packet starts to interact 
with the front pulses. Figure [13] shows typical solutions 
for the different values of 6 by using the physical param- 
eters of W in the second-order model. We performed 
600 different realisations on the initial random condition 
and calculated the separation distances from the first 12 



pulses for 6 — 0.98 and 1.28, and 1200 realisations taking 
into account the first 6 pulses for 5 = 1.82. 

Figure [T4l depicts the results obtained for the first- and 
second-order models at 8 = 0.98 for W and WG at Re = 3 
and 1.59, respectively. Panel (a) shows the histogram of 
the separation distances between pulses obtained from 
the first-order model and panels (b) and (c) from the 
second-order model for W and WG, respectively. The 
inset in each panel shows the corresponding Si and S2 
functions from the theoretical model given by (|30l) . We 
first note that in all the cases the distributions for the 
pulse separation distances are mainly characterized by 
three peaks that correspond in each case to the theo- 
retically predicted distances at which two-pulses bound 
states are formed (see table in Fig. [T4"|) . Furthermore, 
it is interesting to note that the peaks of the distribu- 
tions appear to be more pronounced and sharper for the 
simulations of the second-order model, specially the ones 
observed at short distances, i\ and £2- Also, while the 
locations at which bound-states are formed are approxi- 
mately the same in all cases, the dominant peaks in the 
distribution change as the viscous-dispersion effects are 
increased. In particular, we observe that the distribu- 
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tion for the first-order simulations is mainly dominated 
by the bound states formed at £2, whereas the distri- 
butions for the second-order simulations are dominated 
by the peaks at £\ and £2, particularly in the WG case 
where both peaks have approximately the same proba- 
bilities. We also note a small peak, less than 0.5%, in 
the second-order WG results, reflecting the fact that the 
first observed bound state can indeed be affected by the 
second-order terms. A more clear picture on the effect 
of viscous dispersion on pulse dynamics is expected to 
emerge as S is increased. 

Figure [TS] shows the results for 5 = 1.28, that cor- 
responds to Re = 3.75 and 2 for W and WG, respec- 
tively. The distributions are again characterized by the 
presence of peaks representing the selected distances at 
which bound states are observed. In this case, we ob- 
serve four peaks that are sharper and more pronounced 
in the second-order simulations than we observed before. 
We find that the histogram peaks obtained in the W- 
simulations occur at similar distances for both the first- 
and the second-order models, and are in very good agree- 
ment with the theory (see the table in Fig. [T5"|) . It is in- 
teresting to note, however, that such a similarity between 
both models is no longer observed when the simulations 
are made by using the parameter values for WG. In this 
case, the second-order effects start to become important, 
and the distances at which the bound states are observed 
change accordingly. We find that the first observed peak 
for the second-order model occurs at £\ ~ 17.4, in con- 
trast to the value obtained by using the first-order model, 
£\ ~ 20.9. Also, we note that the dominant peak in 
the histogram for the second-order model is found at 
£2 — 19.9, whereas the peak corresponding to the first- 
order model is found at £3 ~ 31. All of these differences 
can be explained in terms of the decrease of the ampli- 
tude of the capillary ripples preceding the main solitary 
humps that is largely influence by second-order effects (cf. 
Fig. [2]), which in turn, allows the pulses to get closer to 
each other. It is worth mentioning that for £\ w 17.4 the 
assumption of well-separated pulses is not quite satisfied 
and, as expected, the theoretical and numerical results 
are not in good agreement. 

Figure [TBI depicts the results for 5 — 1.82, that corre- 
sponds to Re = 5 and 2.66 for W- and WG-simulations, 
respectively. The W-simulations reveal that there is no 
any clear distance selection among pulses. The distribu- 
tions of the separation distances for both the first- and 
second-order models are broader, as compared to the pre- 
vious cases and without a dominant peak. This seems to 
indicate that the pulses start to feel the effect of some 
underlying chaotic behaviour with the eventual loss of 
self-organisation into bound states, but unlike the spatio- 
temporal chaos with the KS equation, solitary pulses are 
still clearly identifiable. Interestingly, when the simu- 
lations are performed using the values of the parame- 
ters corresponding to WG, the consequence of includ- 
ing the viscous-dispersion effects due to the second-order 
terms becomes crucial in the bound-state formation pro- 





5 = 0.98 


5 


= 1.28 


lst-order 


0.036 




0.034 


2nd-order (W) 


0.038 




0.035 


2nd-order (WG) 


0.039 




0.045 



TABLE I. Density p a of the solitary-wave gas for different 
values of 5 and for the first-order model, and second-order 
model for both W and WG. 



cess. Indeed, the histogram obtained for the second-order 
model is still dominated by the presence of three peaks 
which are in agreement with the theoretically predicted 
bound states. This is a clear evidence that viscous dis- 
persion plays an important role in the self-organisation 
process that brings the system in a state which can be 
described in terms of bound states. 

C. Solitary-wave gas density 

Our numerical observations for 6 = 0.98 and 1.28 have 
shown that there is a clear selection process for which 
the pulses self-organise to form bound states described by 
well-defined distances. From a statistical point of view, it 
is therefore reasonable that in a large domain containing 
many solitary waves interacting with each other, the sys- 
tem can be described in terms of a mean separation dis- 
tance, that gives information about the density of solitary 
waves. In this sense, we can treat the system as a "gas" 
compound of solitary waves with a well-defined density. 
To characterize such a solitary-wave gas, we shall assume 
that for large times and large spatial domains, the dis- 
tribution of the separation distances between pulses is 
mainly dominated by the peaks observed in Figs. 1141 and 
1151 We then define the mean separation distance as 

n 

{£) = J2<*di, (34) 
i=i 

where n is the number of peaks of the histogram, and on 
is an average weight that we approximate as 

&i = „ro ) ( 35 ) 

where Pi is the probability of each peak. By using this 
definition we can estimate the typical mean separation 
length between the pulses and therefore, we can obtain 
an estimate of the density of the solitary-wave gas by 
using: 



The densities obtained for each S are given in Table HI As 
expected, the density of the solitary-wave gas increases 
as the viscous-dispersion effects are taken into account. 

The existence of a well-defined density for the solitary 
wave gas indicates the formation of a "permanent" self- 
organized state. An important feature of this state, the 
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FIG. 16. Histograms of the pulse-separation distances for S — 1.82 in (a) the first-order model (77 = 0), (b) the second-order 
model with the parameter values corresponding to W (77 = 0.015), and (c) the second-order model with the parameter values 
corresponding to WG (77 = 0.052). The table shows the numerical value of the peaks observed in the histograms compared 
to the values predicted by the model given by (|29a|) . plotted in the insets of each panel, where the dashed and solid lines 
corresponds to the & and Si functions, respectively. 



average separation distance between the pulses, is largely 
dependent on viscous-dispersion effects. 



V. CONCLUSIONS 

We have examined both analytically and numerically 
the interactions of two-dimensional solitary pulses in 
falling liquid films. We focused in particular on the for- 
mation of bound states of pulses and how it is affected by 
the second-order (in the long- wave expansion parameter) 
viscous-dispersion effects. To this end, we made use of a 
second-order two-field model derived in Refs. and [Til. 
This is a system of coupled nonlinear partial differential 
equations for the local flow rate and the film thickness. 

Our theoretical investigation of the formation of bound 
states was based on appropriately extending the rigorous 
coherent-structure theory for the gKS equation recently 
developed in Refs. HH, 0, and [HI to the second-order 
two-field model (and thus putting the coherent-structure 
theory for falling films on a rigorous basis). By assum- 
ing that the solution is given by a superposition of N 
pulses and a small overlap (correction) function, i.e. as- 
suming that the pulses are well separated (weak inter- 
action), we were able to write down a dynamic equation 



for the overlap function in the vicinity of each pulse, that 
is described by a linear matrix/differential operator and 
contains forcing terms due to the neighboring pulses. A 
careful and detailed analysis of the spectral properties 
of the linear operator revealed that the relevant eigen- 
function is the translational mode that corresponds to 
the zero eigenvalue. This eigenvalue is not isolated and, 
therefore, belongs to the essential spectrum. The null 
space of the adjoint operator is spanned by a constant 
vector function and another non-constant function that 
has an infinite norm, meaning that zero is not in the point 
spectrum of the adjoint operator. This spectral behavior 
is similar to the one observed for the gKS equatio n 20 ' 21 . 
Projections to the translational mode were made rigorous 
by using formulation in a weighted space. The outcome 
of the projections is a dynamical system for the pulse 
locations. By studying its fixed points, we were able 
to predict the distances at which the bound states are 
formed. 

Numerical experiments of the temporal evolution of a 
superposition of two pulses have been found in very good 
agreement with the theoretical predictions, and in par- 
ticular, the theoretically predicted attractive and repul- 
sive dynamics that gives rise to the formation of bound 
states. We demonstrated that the second-order viscous- 
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dispersion terms are crucial for an accurate description 
of the pulse interactions. These terms affect the ampli- 
tude and frequency of the capillary ripples in front of a 
solitary pulse. So their influence is in fact linear, but in- 
terestingly they can have some crucial consequences on 
the nonlinear dynamics of the film and the wave-selection 
process in the spatio-temporal evolution. After all, soli- 
tary pulses interact through their tails which overlap, i.e. 
the capillary ripples and their amplitude and frequency 
will affect the separation distance between the pulses: for 
example, smaller-amplitude ripples will allow for more 
overlap between the tails of neighboring pulses, thus de- 
creasing their separation distance. This in turn will affect 
the average separation distance between pulses when the 
system reaches its permanent quasi-turbulent regime and 
hence the density of the solitary waves. This also means 
that any model that does not include second-order terms 
should be used with caution and certainly not for an 
accurate description of solitary pulse interaction which 
dominates the spatio-temporal dynamics of the film. 

In addition, we have studied strong interaction be- 
tween two pulses and found that it leads to an oscillatory 
behavior of the separation length as was also numerically 
observed in Ref. i33J from direct numerical simulations of 
full Navier-Stokes equations and wall and free-surface 
boundary conditions. This behavior escapes the descrip- 
tion of our weak-interaction theory. Again, we find that 
viscous dispersion effects are crucial and the strong non- 
linear interaction between pulses can be largely reduced 
by increasing the viscous-dispersion parameter. 

We have also performed numerical simulations in ex- 
tended domains with a localized random initial condition. 
These allowed us to study the interaction between the 
pulses and how they self-organize to form bound states 
in extended domains. Detailed statistical analysis of the 
pulse separation distances revealed that the histograms 
of the separation distances have clear peaks that corre- 
spond to the theoretically predicted bound states. We 
have used different values of the reduced Reynolds num- 
ber 5 and two different sets of physical parameters cor- 
responding to two liquids of different viscosities, namely 
water and a mixture of water and glycerin. In all cases, 
we have observed that the peaks of the histograms corre- 
sponding to the shortest distances are always more pro- 
nounced in the second-order simulations. As expected, 
the differences between the first- and the second-order 
models were found to be more important in the case of 
the higher viscosity liquid. We observed that the mini- 
mum distance at which bound states are formed is always 
shorter in the second-order simulations. It is important 
to note that in the cases of S = 0.98 and 1.28, both mod- 
els have shown that there is always formation of bound 
states, indicating that statistically, the free surface can 
be treated as a solitary-wave gas, characterized by a typ- 
ical constant mean distance between the pulses. 

Of particular interest would be the extension of the 
coherent structures theory developed here to other vis- 
cous flow problems, for example the problem of a thin 



film coating a vertical fiber—— in which case due to 
the Rayleigh-Plateau instability the film breaks up into 
a train of droplike solitary waves. The recent experi- 
ments in Ref. 137J indicate clearly the formation of bound 
states in a region of the parameter space where the 
Rayleigh-Plateau instability competes with viscous dis- 
persion. The qualitative agreement between the exper- 
iments and the coherent-structure theory developed in 
Refs. HO and [2l| is encouraging. Our hope is that quan- 
titative agreement can be achieved by extending the 
present theory to the second-order model for flow down a 
fiber developed in Ref. H3- For that matter, the present 
formalism could be extended to other two-equation sys- 
tems, e.g. coupled KS-type equations used to describe 
synchronization phenomena 38 . 

Finally, as noted in the Introduction, the two- 
dimensional pulses considered here are only observed 
up to a certain, low-to-moderate, value of the Reynolds 
number. At higher Reynolds numbers, two-dimensional 
pulses develop an instability in the transverse direction 
and a transition to a fully develo ped three-dimensional 
regime is observed (e.g. and |H). The coherent- 
structure theory developed here can be viewed as a foun- 
dation first step for the analysis of the interactions be- 
tween three-dimensional pulses in falling films. 
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Appendix A: The adjoint operators 

The adjoint operator C* is defined as 

(u,dv) = {qu,v), (Al) 

where (•, •) denotes the usual inner product in given 
by 

/oo 
u T -vdx. (A2) 
-oo 

After integration by parts, we find that the adjoint oper- 
ator is 

= ( £*3 £*4 ) i (A3) 
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FIG. 17. Spectrum of C* on an infinite domain for 8 — 0.98 
with the physical parameters corresponding W. The solid line 
is the locus of the essential spectrum and the cross represents 
the point spectrum. 
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The spectrum of the adjoint operator is shown in Fig. [17] 
In addition, it is straightforward to see that on a periodic 
domain the zero eigenfunction is a constant: 



(A4) 



so that C**S?\ — 0. As shown in Fig. Effe), on an infi- 
nite domain we have m — > and therefore there is no 
such a function in the null space of the adjoint operator 
that also belongs to Lg. We can therefore conclude that 
zero is not in the point spectrum of C* . The generalized 
zero eigenfunction on a periodic domain, i.e. the func- 
tion satisfying C*^\ = ^f\, is found numerically and its 
components are shown in Fig. [5] 

We can also show that the adjoint operator S* is: 



J* = 



J* 1 
J* 2 



(A5) 
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